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IMAGING OF SCATTERING MEDIA 
USING RELATIVE DETECTOR VALUES 



This invention was made with U.S. Government support under contract number 
5 CA-RO 166 184-02 A, awarded by the National Cancer Institute. The U.S. Government 
has certain rights in the invention. 

This application claims the benefit under 35 U.S.C. §120 of prior U.S. Provisional 
Patent Application Serial Nos. 60/153,769 filed September 14, 1999, entitled 
TOMOGRAPHY IN A SCATTERING MEDIUM, 60/153,926 filed September 14, 1999, 
1 0 entitled DYNAMIC TOMOGRAPHY IN A SCATTERING MEDIUM and 60/154,099 
q filed September 15, 1 999, entitled DYNAMIC TOMOGRAPHY IN A SCATTERING 
m MEDIUM. 

H This application is related to copending application serial number "not vet 

* w assigned", attorney docket number 0887-4147PC 1 , filed on the same date as this 

Hl5 application, entitled "SYSTEM AND METHOD FOR TOMOGRAPHIC IMAGING OF 

Sirs? 

J DYNAMIC PROPERTIES OF A SCATTERING MEDIUM" by inventors R. Barbour 
fU and C. Schmitz and is hereby incorporated by reference (hereinafter the "Barbour 
4147PC1 application") 

This application is also related to copending application serial number "not vet 
20 assigned", attorney docket number 0887-4 149PC1, filed on the same date as this 

application, entitled "METHOD AND SYSTEM FOR IMAGING THE DYNAMICS OF 
A SCATTERING MEDIUM'* by inventor R. Barbour and is hereby incorporated by 
reference (hereinafter the "Barbour 4147PC2 application"). 
Field of the Invention 
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The invention relates generally to imaging in a scattering medium, and more 
particularly, to a method using a novel modification to the perturbation formulation of the 
radiation transport inverse problem to determine relative changes in the absorption and/or 
scattering properties of the medium based on relative changes in measured energy. 
Background of the Invention 

Many techniques and systems have been developed to image the interior structure 
of a turbid medium through the measurement of energy that becomes scattered upon 
being introduced into a medium. Typically, a system for imaging based on scattered 
energy detection includes a source for directing energy into a target medium and a 
plurality of detectors for measuring the intensity of the scattered energy exiting the target 
medium at various locations with respect to the source. Based on the measured intensity 
of the energy exiting the target medium, it is possible to reconstruct an image 
representing the cross-sectional scattering and/or absorption properties of the target. 
Exemplary methods and systems are disclosed in Barbour et al., U.S. Patent No. 
5,137,355, entitled "Method of Imaging a Random Medium," (hereinafter the "Barbour 
'355 patent"), Barbour, U.S. Patent No. 6,081,322, entitled "NIR Clinical Opti-Scan 
System," (hereinafter the "Barbour '322 patent"), the Barbour 4147PC1 application, and 
the Barbour 4147PC2 application. 

Imaging techniques based on the detection of scattered energy are capable of 
measuring the internal absorption, scattering and other properties of a medium using 
sources whose penetrating energy is highly scattered by the medium. Accordingly, these 
techniques permit the use of wavelengths and types of energy not suitable for familiar 
transmission imaging techniques. Thus they have great potential for detecting properties 
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of media that are not accessible to traditional energy sources used for transmission 
imaging techniques. For example, one flourishing application of imaging in scattering 
media is in the field of optical tomography. Optical tomography permits the use of near 
infrared energy as an imaging source. Near infrared energy is highly scattered by human 
tissue and is therefore an unsuitable source for transmission imaging in human tissue. 
However, these properties make it a superior imaging source for scattering imaging 
techniques. The ability to use near infrared energy as an imaging source is of particular 
interest in clinical medicine because it is exceptionally responsive to blood volume and 
blood oxygenation levels, thus having great potential for detecting cardiovascular disease, 
tumors and other disease states. 

A common approach for the reconstruction of an image of the cross-sectional 
properties of a scattering medium is to solve a perturbation equation based on the 
radiation transport equation. The radiation transport equation is a mathematical 
expression describing the propagation of energy through a scattering medium. The 
perturbation formulation relates the difference between coefficient values of the true 
target and a specified reference medium, weighted by a proportionality coefficient whose 
value depends on, among other things, the source/detector configuration and the optical 
properties of the medium. In practice, tomographic measurements consider some array of 
measurement data, thus forming a system of linear equations having the form 



where 6u is the vector of differences between a set of measured light intensities (u) and those 
predicted for a selected reference medium (u r ), W r is the Jacobian operator, and Sx is the 
position-dependent difference between one or more optical properties of the target and reference 
media (i.e., a change in absorption coefficient Sfi* a change in the reduced scattering coefficient, 



u-u r = Su = W r <£x, 



(1) 
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/f* or, in the diffusion approximation, the diffusion coefficient SD, where Z)= l/[3(p a + //',)]). 
The operator, referred to as the weight matrix, has coefficient values that physically represent the 
fractional change in light intensity at the surface caused by an incremental change in the optical 
properties at a specified point in the medium. Mathematically this is represented by the partial 
5 differential operator dujdxj, where i is related to the i* source/detector pairs at the surface of the 
medium, andj to the 7 th pixel or element in the medium. 

Although the perturbation equation in Eq. (1) can be solved using any of a number of 
available inversion schemes, practical experience has shown that the accuracy and reliability of 
the results obtained can vary greatly due to uncertainties and errors associated with the quality of 
10 the measurement data, inaccuracies in the physical model describing light propagation in tissue, 
Q specification of an insufficiently accurate reference state, the existence of an inherently 
CO underdetermined state caused by insufficiently dense measurement sets, weak spatial gradients in 
M* the weight function, and so forth. 

Ill In Practice, a matter of considerable concern is the accuracy with which the reference 

fc§ medium can be chosen. An accurate reference is one that closely matches the external geometry 
jj. of the target medium, has the same size, nearly the same internal composition, and for which the 
Q ^cations of the measuring probes and their efficiency coincide well with those used in the actual 

PS 5 

I y 

measurements. While such conditions may be easily met in numerical and perhaps laboratory 
phantom studies, they represent a much greater challenge in the case of tissue studies. 

20 Confounding factors include the plasticity of tissue (it deforms upon probe contact), its mainly 
arbitrary external geometry and internal composition and the considerable uncertainty stemming 
from the expected variable coupling efficiency of light at the tissue surface. The influence of 
these uncertainties can be appreciated when it is recognized that the input data vector for the 
standard perturbation formulation (£*., Eq. (1)) is actually the difference between a measured and 

25 a computed quantity. This vector contains information regarding the subsurface properties of the 
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target medium that, m principle, can be extracted provided an accurate reference medium is 
available. 

In practice, however, there are two significant concerns that are frequently encountered in 
experimental studies and are not easily resolvable especially in the case of tissue studies. One 
5 concern is the expected variable coupling efficiency of light entering and exiting tissue. 
Nonuniformity in the tissue surface, the presence of hair or other blemishes, its variable 
deformation upon contact with optical fibers, the expected variable reactivity of the vasculature in 
the vicinity of the measuring probe all serve to limit the ability to accurately determine the in- 
coupling and out-coupling efficiencies of the penetrating energy. Consideration of this issue is 
10 critical as variations in the coupling efficiency will be interpreted by the reconstruction methods 
3 as variations in properties of the target medium and can introduce gross distortions in the 

VS. 

q recovered images. In principle, the noted concern can be minimized by adopting absolute 
calibration schemes, however, in practice the variability in tissue surface qualities will limit 
reliability and stability of these efforts. 
ff|5 A second concern stems from the underlying physics of energy transport in highly 

£T scattering media. One effect of scattering is to greatly increase the pathlength of the propagating 
g energy. Small changes in the estimated absorption or scattering properties of the medium can, 
depending on the distance separating the source and detector, greatly influence the density of 
emerging energy. This consideration has important implications regarding the required accuracy 
20 by which the reference medium must be specified. In the context of perturbation formulations, 
the reference medium serves to provide estimates of the predicted energy density as well as to 
provide the needed weight functions that serve as the imaging operators. The difficulty is that the 
computed reference intensity values are extremely dependent on the optical coefficient values of 
the reference medium. Significantly, this dependence is a nonlinear function of the distance 
25 between source and detector. It follows that a small change in the optical properties of the 
reference medium may influence the value of the computed intensity differences (<*) by a 
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relative amount that may be significantly different for each source/detector pair, thereby altering 
the information content of the data vectors. This can lead to the recovery of grossly corrupted 
images. Whereas, in principle, such effects may be overcome by use of recursive solutions to the 
perturbation equation (i.e., Newton-type updates), in practice this can require extensive 
5 computational efforts, especially in the case of 3D solutions. Moreover, it is well appreciated that 
such efforts to improve on first order solutions to the perturbation equation (e.g., Born or Rytov 
solutions), can fail if the initial estimate chosen for the reference medium is insufficiently 
accurate. 

One alternative to devising absolute calibration schemes is to devise 
10 methodologies whose solutions are intrinsically less sensitive, or better still, do not 
M require such information, but nevertheless are capable of providing accurate descriptions 
S ° f certain features of highly scattering media. While a range of empirical methodologies 
h can be devised » ^ is desirable that they be broadly extendable without requiring undue 

physical approximations, since these are generally incompatible with model-based 
hi 5 methods. 

jF An approach previously adopted is to directly relate relative detector readings, 

lU obtained from comparison of detector values derived from two different target media 
(usually media with and without the included object), to the weight matrix computed 
based on a previously assigned reference medium. R. L. Barbour, H. Graber, R. Aronson, 
20 and J. Lubowsky, "Model for 3-D optical imaging of tissue," Int. Geosci. and Remote 
Sensing Symp., (IGARSS), 2, 1395-1399 (1990). While capable of producing good 
quality images of internal structure of a target medium, the method proved to have 
limited utility as it did not produce solutions having physical units, thereby rendering 
specific interpretation difficult, as well as limiting efforts to compute recursive solutions. 
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For the forgoing reasons, there is a need for image reconstruction techniques 
based on the detection of scattered energy that (1) do not require absolute calibration of, 
and absolute measurements by, the detectors and other elements of the apparatus, (2) 
make the standard perturbation equation less susceptible to variations between boundary 
5 conditions and properties of the reference medium and the target medium, and (3) 
produce solutions having physical units. 

SUMMARY 

The present invention satisfies these needs by providing a method for generating 
an image of a scattering medium using normalized relative measured intensity values a 
Q0 perturbation formulation based on the radiation transport equation, 
w It is an object of the present invention to provide a method for imaging the 

% properties of a scattering target medium using a modified perturbation equation. The 

ro , , 

method comprises generating a first data vector of measured data from a target and a 
y second vector of measured data from a target, normalizing the first and second vectors of 
J; 5 measured data and solving a modified perturbation equation for the unknown optical 
fy properties of a target medium. The first and second vectors of measured data are 
measures of energy emerging from the target. 

It is a further aspect of this invention to obtain the first and second sets of 
measured data from the same target, wherein the first set of measured data is a set of data 
20 measured at an instant in time and the second set of measured data is a time average 
mean of a plurality of first sets of measured data. 

It is yet a further aspect of this invention to obtain the first and second sets of 
measured data from two different targets. 
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In another aspect of this invention the modified perturbation equation is a 
modified Rytov approximation. 

In another aspect of this invention the modified perturbation equation is a modified Born 
approximation. 

BRIEF DESCRIPTION OF THE FIGURES 

For a better understanding of the invention, together with the various features and 
advantages thereof, reference should be made to the following detailed description of the 
preferred embodiment and to the accompanying drawings, wherein: 

FIG. 1 is a schematic illustration of an exemplary imaging system; 

FIG. 2A is a cross-sectional image of the absorption coefficients of the target; 

FIG. 2B is a cross-sectional image of the diffusion coefficients of the target; 

FIG. 3 is a table illustrating a summary of the test cases explored; 

FIG. 4A is a series of reconstructed cross-sectional absorption profile images of 
the target obtained in test case 1 of FIG. 3; 

FIG. 4B is a series of reconstructed cross-sectional diffusion profile images of the 
target obtained in test case 1 of FIG. 3; 

FIG. 5A is a series of reconstructed cross-sectional absorption profile images of 
the target obtained in test case 2 of FIG. 3; 

FIG. 5B is a series of reconstructed cross-sectional diffusion profile images of the 
target obtained in test case 2 of FIG. 3; 

FIG. 6 A is a series of reconstructed cross-sectional absorption profile images of 
the target obtained in test case 3 of FIG. 3; 
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FIG. 6B is a series of reconstructed cross-sectional diffusion profile images of the 
target obtained in test case 3 of FIG, 3; 

FIG. 7 A is a series of reconstructed cross-sectional absorption profile images of 
the target obtained in test case 5 of FIG. 3; 

FIG. 7B is a series of reconstructed cross-sectional diffusion profile images of the 
target obtained in test case 5 of FIG. 3 ; 

FIG. 8 A is a series of reconstructed cross-sectional absorption profile images of 
the target obtained in test case 6 of FIG. 3; 

FIG. 8B is a series of reconstructed cross-sectional diffusion profile images of the 
target obtained in test case 6 of FIG. 3; 

FIG. 9A is a series of reconstructed cross- sectional absorption profile images of 
the target obtained in test case 7 of FIG. 3; 

FIG. 9B is a series of reconstructed cross-sectional diffusion profile images of the 
target obtained in test case 7 of FIG. 3; 

FIG. 10 is a table listing constant calibration errors corresponding to each image 
of FIGS. 9A and 9B; 

FIG. 1 1 is a table of data corresponding to the error, resolution and contrast for 
the reconstructed images shown in FIGS. 4 A and B; 

FIG. 12 is a table of data corresponding to the error, resolution and contrast for 
the reconstructed images shown in FIGS. 5 A and B; 

FIG. 13A is a series of reconstructed cross-sectional absorption profile images of 
the target shown in FIG. 15 with varying concentrations of hemoglobin; 
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FIG. 13B is a series of reconstructed cross-sectional diffusion profile images of 
the target shown in FIG. 15 with varying concentrations of hemoglobin; 

FIG. 14A is a series of reconstructed cross-sectional absorption profile images of 
the target shown in FIG. 15 with varying reference medium properties; 

FIG. 14B is a series of reconstructed cross-sectional diffusion profile images of 
the target shown in FIG. 15 with varying reference medium properties; 

FIG. 15 is a schematic illustration of a phantom study; 

FIG. 16 A is a graph plotting the amplitude of normalized intensities used in the 
modified perturbation formulation corresponding to the reconstructed cross-sectional 
images shown in row 3 of FIG. 4; 

FIG. 16B is a graph plotting the amplitude of normalized intensities used in the 
modified perturbation formulation corresponding to the reconstructed cross-sectional 
images shown in column 3 of FIG. 4; 

FIG. 17A is a graph plotting the amplitude of the relative intensities used in the 
standard perturbation formulation corresponding to the reconstructed cross-sectional 
images shown in row 3 of FIG. 5; 

FIG. 17B is a graph plotting the amplitude of the relative intensities used in the 
standard perturbation formulation corresponding to the reconstructed cross-sectional 
images shown in column 3 of FIG. 5; 

FIG. 18 is a graph plotting the amplitude of the frequency spectrum of the Fourier 
transforms for data vectors computed using equation (1) and equation (3); and 

FIG. 19 is a table illustrating the ratio of average contrasts of reconstructed 
absorption and diffusion coefficients shown in FIGS. 4A and B respectively. 
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DETAILED DESCRIPTION OF THE INVENTION 

Imaging in a scattering medium relates to the methods and techniques of image 
the internal properties of a medium based on the detection of scattered energy. 
The typical process for imaging of a scattering medium comprises: (1) selecting a 
reference medium having known boundary conditions and optical properties which are 
substantially similar to those of the intended target; (2) determining a weight matrix and 
an intensity of emerging energy exiting the reference medium at each of a plurality of 
source points for each of a plurality of detectors located around the reference medium 
boundary, the determination being made by either actual measurements or solution of the 
radiation transport equation; (3) measuring actual emerging energy intensities for 
corresponding source and detector points on a target medium; and (4) solving the 
perturbation equation for the optical properties of the target based on the measured 
intensities of energy emerging from the target. 

The present invention describes an improved methodology for imaging of a 
scattering medium using a modified fomi of the standard perturbation equation. The 
inventive modification of the standard perturbation equation is capable of (1) reducing 
the sensitivity of the perturbation equation to differences between the reference medium 
and target medium, (2) producing solutions to the perturbation equation having physical 
units, and (3) reducing the effect of variable detector efficiencies without the need for 
absolute calibration, while at the same time preserving the ability to compute recursive 
solutions. Compared to the standard perturbation approach, the described invention 
provides remarkable improvement in the quality of image reconstruction. 
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While the method of the present invention is applicable to known static imaging 
techniques it is instrumental in the realization of practical dynamic imaging of highly 
scattering media. There are three principal elements to practical dynamic imaging. The 
first element is the development of a fast, parallel, multi-channel acquisition system that 
employs geometrically adaptive measuring heads. This system is described briefly below 
and in further detail in the copending Barbour 4147PC1 application. The second element 
is to evaluate the acquired tomographic data using the modified perturbation methods of 
the present invention. The third element is to collect a time series of data and subject 
either the time series of data or a time series of reconstructed images from the data to 
analysis using various linear and nonlinear time-series analysis methods to extract 
dynamic information and isolated dynamic information. These methods are described in 
detail in the copending Barbour 4147PC2 application. 

Some of the methods, systems and experimental results described below focus on 
optical tomography of human tissue using wavelengths in the near infrared region for the 
imaging source. However, as disclosed generally herein, it will be appreciated to those 
skilled in the art that the invention is applicable to the use of essentially any energy 
source (e.g., electromagnetic, acoustic, and the like) on any scattering medium (e.g., body 
tissues, oceans, foggy atmospheres, earth strata, industrial materials) so long as diffusive 
type mechanisms are the principal means for energy transport through the medium. 
System 

Numerous imaging systems such as those disclosed in the previously mentioned 
the Barbour '355 patent, the Barbour '322 patent and the Barbour 4147PC1 application 
have been developed for use in imaging of a scattering medium. A schematic illustration 
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of an exemplary system is shown in FIG. 1, This system includes a computer 102, 
sources 104, 106, a source demultiplexer 108, an imaging head 1 10, detectors 1 12 and a 
data acquisition board 114. 

A target 1 16 placed in the imaging head 1 10 is exposed to optical energy from 
sources 104, 106. The optical energy originating from sources 104, 106, is combined by 
beam splitter 1 18 and is delivered to source demultiplexer 108. The source demultiplexer 
108 is controlled by computer 102 to direct the optical energy to source fibers 120 
sequentially. 

Each source fiber 120 carries the optical energy from the demultiplexer 108 to the 
imaging head 1 10, where the optical energy is directed into the target 116. The imaging 
head 1 10 contains a plurality of source fibers 120 and detector fibers 122 for transmitting 
and receiving light energy, respectively. Each source fiber 120 forms a source-detector 
pair with each detector fiber 122 in the imaging head 1 10 to create a plurality of source- 
detector pairs. The optical energy entering the target 1 16 at one location is scattered and 
may emerge at any location around the target 1 16. The emerging optical energy is 
collected by detector fibers 122 mounted in the imaging head 1 10. 

The detector fibers 122 carry the emerging energy to detectors 1 12, such as 
photodiodes or a CCD array, that measure the intensity of the optical energy and deliver a 
corresponding signal to a data acquisition board 114. The data acquisition board 1 14, in 
turn, delivers the data to computer 102. 

This imaging process is repeated so as to deliver optical energy to each of the 
source fibers sequentially, a measurement being obtained for detected emerging energy at 
each detector for each emitting source fiber. This process may continue over a period of 
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time with the computer 102 storing the data for reconstruction of one or more images. 

Additionally, the system may include two or more imaging heads for comparing one 

target to another. The computer 102 reconstructs an image representative of the internal 

optical properties of the target by solving a perturbation equation. It will be appreciated 
5 by those skilled in the art that more than one computer can be used to increase data 

handling and image processing speeds. 

The Standard Perturbation Formulation 

As discussed above, reconstruction of a cross section image of the absorption 

and/or scattering properties of the target medium is based on the solution of a 
WO perturbation formulation of the radiation transport equation. The perturbation method 
2 assumes that the composition of the unknown target medium deviates only by a small 
fz amount from a known reference medium. This reduces a highly non-linear problem to 
Jfi one that is linear with respect to the difference in absorption and scattering properties 
q between the target medium under investigation and the reference medium. The resulting 
Ml 5 standard perturbation equation has the following forms: 



In equation (1), 5u is a vector of source-detector pair intensity differences between the 
measured target medium and the known reference medium (i.e., 5u = u - u r ). W is the 
weight matrix describing the influence that each volume element ("voxel") of the 
20 reference medium has on energy traveling from each source to each detector, for all 
source-detector pairs. The volume elements are formed by dividing a slice of the 
reference medium into an imaginary grid of contiguous, non-overlapping pieces. 
Physically, the weight matrix contains the first order partial derivatives of the detector 



u-u r = Su = W r <5x, 



(i) 
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responses with respect to the optical coefficients of each volume element of the reference 
medium. 5x is the vector of differences between the known optical properties (e.g., 
absorption and scattering coefficients) of each volume element of the reference medium 
and the corresponding unknown optical properties of each volume element of the target 
medium. 

This standard perturbation equation assumes (1) use of absolute detector 
measurements for u, and (2) that the boundary conditions and optical properties of the 
reference do not vary significantly from the target. Both of these factors are problematic 
in practice. 

The Modified Perturbation Formulation 

The present invention modifies the standard perturbation equation by replacing 5u 
with a proportionate relative difference between two measured values multiplied by a 
reference term of the required units as set forth in the equation (2) below: 



where i is the source/detector pair index. In equation (2), I r is the computed detector 
reading corresponding to a source-detector pair of a selected reference medium, and I and 
Io represent two data measurements for a corresponding source-detector pair on one or 
more targets (e.g., background vs. target, or time-averaged mean vs. a specific time point, 
etc.). The resultant term 8I r therefore represents a relative difference between two sets of 
measured data that is then mapped onto a reference medium. Careful examination 
reveals that this modification has important attributes that limit the effects of modeling 



{*,). - 



(Io),- 



(U 



(2) 
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errors and minimize ill-conditioning of the inverse problem while retaining the correct 
units in the solution. 

The corresponding perturbation equation using this modified term is: 

w r .^x = ^l r (3) 

In equation (3) W r and 5x are the same as W r and 5u in equation (1). Referring to 
equations (2) and (3), it can be seen that in the limit, when I r equals to I 0 , this equation 
reduces to the standard perturbation formulation shown in equation (1). This formulation 
represents the Born approximation formulation of the modified perturbation equation. A 
similar expression may be written for the Rytov approximation in the following form: 

(wA. 

(w;) ^7u ; (4) 

The inventive operation accomplished by equation (2) is to preserve the information 
content of a measured proportionate relative data vectors obtained from the target 
medium and to yield data vectors having the correct physical units. Apart from 
simplifying measurement requirements, the method represented by equations (3) and (4) 
also reduces susceptibility of the perturbation equation to boundary and optical property 
variation between the target and the reference medium. Additionally, iterative solutions 
of equations (3) and (4) can be easily implemented. In principle, the techniques used in 
the modified perturbation equation, referred to as the normalized difference* method 
(NDM), can be used to evaluate any measured proportionate relative quantity. 
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Experimental Validation 

The following discussion presents results validating the methods and advantages 
of the present inventions. These examples are presented merely as an illustration of the 
benefits of applying the NDM approach for the analysis of relative measures from highly 
scattering media. 

Forward Model and Data Acquisition Geometry 

For any intended target the perturbation equation is generated for a reference 
medium having boundary conditions and optical properties substantially similar to the 
target. The perturbation equation models the energy propagation, e.g. light, in the 
reference medium as a diffusion process. For a domain Q having a boundary dCl, this is 
represented by the expression: 



where u(r) is the photon density at position r, r s is the position of a DC point source, and 
D(r) and fJL Q (v) are the position-dependent diffusion and absorption coefficients, 
respectively. The diffusion coefficient is defined as: 



where ft r s [v] is the reduced scattering coefficient. The photon density values at the 

detectors, i.e., the calculated energy intensity emerging from the reference medium at 

each detector, were computed by applying Dirichlet boundary conditions on an 

extrapolated boundary. Depending on the target medium to be explored, sources and 

* 

detectors for the reference are positioned 1 to 2 transport mean free pathlengths within 
the boundary of the reference medium. 



V.[D(r)V«(r)]- A „(r)«(r) = ^(r-r J ), reQ 



(5) 




(6) 
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Solutions to the diffusion equation may be computed by any known means, such 
as by the KASKADE adaptive finite element method. R. Beck, R. Erdmann and R. 
Roitzsch, "Kaskade 3.0 - An object-oriented adaptive finite element code," Technical 
report TR 95-4, Konrad-Zuse-Zentrum fur Informationstechnik, Berlin (1995). This is a 
publicly available code suitable for the solution of partial differential equations in one, 
two or three dimensions using adaptive finite element techniques. The code can be 
readily modified to permit solutions to the diffusion equation using a point source. Mesh 
generation may be by any known method, such as the Delaunay tessellation algorithm 
originally proposed by Watson. D. F. Watson, "Computing the n-dimensional Delaunay 
tessellation with applications to Voronoi polytopes", Computer Journal, 24, 167-172 
(1981). 

The perturbation equation is specific to the boundary conditions and optical 
properties of the reference medium, including the orientation of the source-detector pairs 
in relation to one another and the reference medium. These conditions and properties are 
preferably nearly identical to the target. For example, in the experiments discussed 
below, the perturbation equation was generated based on an imaging system having six 
sources and eighteen detectors per source (108 source-detector pairs) with the sources 
equally spaced at 60 degree intervals around the boundary of the medium and the 
detectors equally spaced at 20 degree intervals. 
Inverse Algorithm 

As described above, in the present invention relative intensities are measured for 
all source-detector pairs using any known imaging system. Image recovery is then 
achieved using known methods, such as conjugate gradient descent (CGD), or 
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simultaneous algebraic reconstruction techniques (SART), to solve the modified 
perturbation equation for the absorption and scattering properties of the target. J, Chang, 
H. L. Graber, R, L. Barbour and R. Aronson, "Recovery of optical cross-section 
perturbations in dense-scattering media by transport-theory-based imaging operators and 
steady-state simulate data", Appl. Opt. 35, 3963-3978, (1996) (the disclosure of which is 
incorporated herein by reference). For example, the experimental results discussed below 
were achieved using a CGD solver with and without matrix rescaling. In addition, a 
weight matrix rescaling (WMR) technique may be used to improve the ill-conditioning of 
the weight matrix. The effect of rescaling the weight matrix is to make it more uniform. 
Two rescaling criteria can be applied for this purpose: (1) rescaling the maximum of each 
column to 1 ; or (2) rescaling the average of each column to 1 . In the experimental results 
below, when WMR was used, criterion 1 was applied for image recovery. 
The solution to the modified perturbation equation provides a relative measure of the 
difference between the cross-sectional optical properties of a target during the first and 
second measurements I and I<>. The values from this solution are used to generate cross- 
sectional images representative of the target's internal optical properties. 
Test Cases Explored 

The following discussion presents results obtained for seven test cases comparing 
image reconstruction using the known standard perturbation formulation with the 
modified perturbation formulation of the present invention. These examples are 
presented merely as an illustration of the benefits of the modified perturbation method of 
the present invention. 
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The reconstruction results presented in the test case are limited to solution of the 
first order Born approximation. The coefficient values for absorption and diffusion 
coefficients were computed simultaneously. For each case tested, measures of error, 
contrast accuracy and resolution were also computed. These are defined as follows: 



where aj and xj are the actual and reconstructed values of the optical coefficient, and M is 
the number of volume elements used for reconstruction. 

ii. Image contrast accuracy was determined by computing the mean value of 
the recovered perturbation coefficient along the transect bisecting the two objects. 

iii. Resolution was measured by computing the mean value of the full-width 
half-maximum of the two reconstructed objects along the transect bisecting the 
inclusions. 

FIGS. 2 A and 2B show the cross-sectional geometry and absorption (FIG. 2 A) 
and diffusion (FIG. 2B) coefficient profiles of the target medium explored. The target 
medium is 8cm in diameter and has two included objects each 1 cm in diameter and 
separated by 3 cm symmetrically about the center of a homogeneous background 
medium. Optical properties of the background and included objects are 0.04 and 0.02 
cm -1 for fii a (absorption coefficient), and 10 and 5 cm" 1 for p s (scattering coefficient), 
respectively. 



Image error is the relative root mean square error (RMSE) 
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The table of FIG. 3 lists the various test cases explored. The symbol "V" 
indicates that the parameter was varied, "C" indicates that the parameter was held 
constant, and *7" indicates that the parameter was not considered. The test cases allowed 
at least partial isolation of the effect of variations in each of the input parameters on the 
resultant image for different reconstruction schemes and perturbation formulations. This 
testing permitted exploration of the dependence of ill-conditioning on the different input 
parameters that influence the accuracy and stability of the image reconstruction. 

Test cases 1 and 2 examined the general case where the reference medium is 
based only on an estimate of the background optical properties of the target medium. The 
estimated properties were varied over a broad range, ranging from 0.0 cm" 1 to 0.3 cm" 1 in 
Ha*, and from 3 cm* 1 to 30 cm" 1 in pi s - For purposes of comparison, test cases 3 through 7 
explored the dependence of image quality on the varied parameters using the standard 
perturbation formulation. 

Test cases 3 and 4 mainly mirror conditions explored in cases 1 and 2 with the 
exceptions that the standard perturbation formulation was evaluated, and a narrower 
range of coefficient values was considered for the reference medium. Here the general 
case is also considered where only an estimate of the background optical properties of 
the target medium is available. The range of values for the optical properties explored 
were from 0.02 cm" 1 to 0.08 cm" 1 in and from 5 cm" 1 to 15 cm" 1 in 

Cases 5 and 6 consider the special situation where prior knowledge of the 
background properties of the target medium is known. The parameters varied were W r 
and I r , referred to as Wb and lb, respectively. The range of optical properties varied for 
test 5 is same as in case 1. For case 6, the range of optical properties varied is the same 
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as in case 3. Test case 7 explores the effect of a constant calibration error in 
measurement, and assumes prior knowledge of the background properties of the target 
medium. 
Results 

Data shown in FIGS. 4 through 9 illustrates the influence that the varied 
parameters listed in FIG. 3 have on the reconstruction results derived from a first-order 
Born approximation using the standard and modified perturbation formulations. The 
results presented are listed in a matrix format. The value of the absorption and scattering 
coefficient for the reference medium is fixed for each row and column, respectively. 
Varied is the value of these parameters along the orthogonal direction. Shown in the 
figures are the reconstruction profiles for all test cases explored except case 4, whose 
findings are reported in the text. 
Qualitative Analysis 

Data in FIGS. 4 and 5 show the quality of reconstructed images obtained using 
equation (3). FIGS. 4A and 5 A illustrate the computed absorption maps, while FIGS. 4B 
and 5B show the computed diffusion maps. Compared to the original profiles shown in 
FIG. 2, these findings illustrate that qualitatively accurate results, revealing the two- 
object structure, are obtained over a broad range of values for the selected reference 
medium. Artifact dominated results are limited to cases in the lower right corner of the 
matrix. These correspond to those reference media having absorption and scattering 
coefficient values significantly greater than the background of the target medium. 
Quantitative analysis of these results is presented in the next section. Comparison of 
images shown in FIGS. 4 and 5 reveals that the matrix rescaling method is capable of 



22 



WO 01/20546 PCT/USOO/25156 

providing a higher resolution image, though over a reduced range of values for the 
reference medium. For example, comparison of results reveals that the matrix rescaling 
method yields only artifacts in the absorption map for non-absorbing reference media, 
while under the same conditions the diffusion coefficient maps reveal two completely 
resolved objects. In the absence of matrix rescaling, both coefficient maps reveal the 
presence of the included objects, though with reduced edge resolution and more artifacts 
in the diffusion map. The added improvement using matrix rescaling, however, is 
achieved under the limiting conditions of a range constraint (positive for D and negative 
for Ha)- Overall, the range of values for the reference medium's optical properties for 
which qualitatively accurate maps are obtained for both inverse methods is much greater 
than previously reported. S. R. Arridge, M. Schweiger, "Sensitivity to prior knowledge 
in optical tomographic reconstruction", in Proc. Optical tomography, photon migration 
and spectroscopy of tissue and model media: Theory, human studies, and 
instrumentation, SPIE, 2389, 378-388, (1995) (the disclosure of which is incorporated 
herein by reference). 

Shown in FIG. 6 are results from case 3 evaluated using equation (1). Compared 
to results shown in FIGS. 4 and 5, a more limited range of values of optical properties for 
the reference medium were examined, since outside of this range, only artifact was 
recovered. Even within the explored range, significant instability was observed for 
relatively small variations in the reference medium. This sensitivity indicates a state of 
ill-conditioning that is alleviated using the modified perturbation formulation (equation 
3). Not shown are results from case 4 using the matrix rescaling method. In case 4, even 
greater instability was observed than in the case using CGD only. 
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The parameters varied in the above figures include both the computed reference 
intensity and the weight matrix. This is the general case where both quantities can only 
be estimated and are computed from a specified reference medium. 

Results shown in FIGS. 7 through 9 explore the special cases where errors occur 
only in one parameter. Data in FIG. 7 illustrates the influence of errors in the estimated 
weight matrix. Assumed is prior knowledge of the reference detector intensity values for 
the background medium. In practice this would correspond to situations where a 
measurement was made in the presence and absence of an included object. Inspection 
reveals that qualitatively accurate results are obtained over a significantly broader range 
of reference values than those presented in FIG. 7. This finding suggests that the 
principal origin of the ill-conditioning of the inverse problem is associated with errors in 
the estimated reference intensity. 

Results shown in FIG. 8 explore this possibility directly. In this situation, we 
assume the unlikely case where accurate prior knowledge of the weight matrix for the 
target medium is available. Comparison of the results in FIG. 8 to the results in FIG. 6 
illustrates some improvement in the range of reference media yielding qualitatively 
accurate results. However, this range is small compared to the situations tested in FIG. 8 
for the standard perturbation formulation and for the modified perturbation formulation in 
FIGS. 4 and 5. 

Finally, in FIG. 9, we test the case where accurate prior knowledge for the weight 
matrix and reference detector values are available and a constant error of measurement is 
introduced. The error added to the measured detector reading, I, varies from -50 to 
900%. Variations in constant calibration error corresponding to FIG. 8 are shown in FIG. 
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9. Under these conditions, the results show that a constant calibration error does not 
significantly affect the qualitative accuracy of the computed images. It is worth noting 
that errors of this type will not exist using the modified perturbation formulation. 
Quantitative Analysis 

A quantitative analysis of the image data was made by computing measures of 
error, resolution and contrast. Results shown in FIGS. 1 1 and 12 are the corresponding 
values computed from data shown in FIGS. 4 and 5. The format of the data is the same 
as in FIG. 4 (i.e.,, data in the rows are derived for a fixed value of absorption, while 
column data is derived for a fixed value of scattering). 

Inspection of FIGS. 1 1 and 12 reveals, not surprisingly, that the lowest image 
RMSE is achieved when the selected reference medium matches the background optical 
properties of the target medium. Consistent with this finding is that nearly equivalent 
error values were obtained for those maps in which one of the fixed coefficient values for 
the reference medium matched that of the background. Interestingly, in the absence of 
WMR, whereas these conditions produced the lowest total error values for the image 
map, improved accuracy of object contrast was obtained using reference media generally 
having reduced scattering and increased absorption values compared to those of the 
background. An exception to this was the case where the background absorption level in 
the reference medium was reduced while scattering value matched the object. With 
WMR, the best accuracy for object coefficient values was achieved when the reference 
medium matched the coefficient values of the background medium. This is not 
unexpected given the imposed constraint. Overall there is evidence of a trade-off 
between artifact levels and accuracy in object contrast. Potentially significant is the 
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observation of several instances where considerably enhanced object contrast is seen 
without undue degradation of image quality. This is observed with either inverse 
algorithm, though the trends are somewhat different. In the absence of WMR, enhanced 
contrast for both absorption and diffusion maps are seen using reference media whose 
absorption and scattering coefficient values are lower than the background. Enhanced 
contrast is also seen in the case of a non-absorbing reference medium, although increased 
artifacts are present. With WMR, improved object contrast is also seen with reference 
media having reduced scattering, but the trend in contrast enhancement is opposite for the 
different coefficients and depends strongly on the value of the absorption coefficient. 
Increasing the absorption coefficient value for the reference medium increases the object 
contrast value for absorption while reducing the contrast for the diffusion coefficient. 
Inspection of results reported for edge resolution reveal an under and over estimate of 
object diameter for images computed with and without the WMR, respectively. Whereas 
it is often best to avoid errors in edge resolution, the resolution of the images obtained 
using the WMR method is striking. 

Experimental Validation of Modified Perturbation (NDM) Formulation 

Experimental verification demonstrating that the modified perturbation 
formulation is capable of resolving internal structure of a dense scattering medium is 
given in FIGS. 13 and 14. Tomographic measurements were performed at 780 nm using 
the IRIS imaging system previously described. R.L. Barbour, R. Andronica, Q. Sha, HI. 
Graber, and I. Soller, "Development and evaluation of the IRIS-OPTIsoanner, a general- 
purpose optical tomographic imaging system." In Advance in Optical Imaging and 
Photon Migration, J.G. Fujimoto et al, ed., Vol 21 of OSA Proceeding Series, pp. 251- 
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255, (1998) (the disclosures of which are incorporated herein by reference). The target 
medium was a latex laboratory glove filled with 2% (v/v) Intralipid suspended from a 
holder in a pendant position. Added to the glove were two 1 cm diameter plastic tubes 
filled with varying concentrations of hemoglobin (Hb) in the amount of 5fxm, lO^m, 
20|im, and 40(im. A cross section of the phantom set-up is shown in FIG. 1 5. The pass- 
through diameter of the IRIS imaging head was closed until gentle contact with the glove 
was achieved. The diameter of the glove in the measurement plane was 6.7 cm. Above 
and below this plane, the glove assumed an arbitrary geometry. Tomographic 
measurements were performed using the same measurement geometry described for the 
numerical studies. Optical measurements were performed in the presence and absence of 
the included objects from which the relative intensity values were derived. The resultant 
data vectors were then evaluated by equation (3) using a regularized CGD method 
without weight matrix rescaling. The coefficient values for the reference medium used 
were varied from 0.01 to 0.04 cm" 1 in fi a and 10 to 20 cm" 1 in // 5 . 

Reconstructed /i a and D maps for each experiment with different concentrations of 
hemoglobin in the two tubes using a specific reference medium are shown in FIG. 13. 
Inspection reveals two well-resolved objects whose contrast in both coefficients increases 
with increasing concentration of added absorber. Quantitatively, the dependence of 
image contrast on absorber concentration is less than expected (bju a recovered = .005 vs. 
bfia actual = .015 cm" 1 ) indicating perhaps either a self-shielding effect, limits of a first 
order Born solution or both. Results shown in FIG. 14 demonstrate that images of similar 
quality are derivable over a range of reference media, a result consistent with the above 
described numerical studies. 
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Discussion and Conclusions 

The present invention describes and evaluates a new formulation for the inverse 
problem for imaging in highly scattering media. Motivating this development has been 
an appreciation of the expected limits imposed by practical measurements, especially as it 
relates to the dependence of image accuracy on instrument calibration and the ability to 
specify an accurate reference medium. This concern arises because many of the 
anticipated clinical applications will require some level of accuracy in the computed 
coefficient values. An accurate solution will require an explicit accounting of various 
factors intrinsic to the detector (e.g., quantum efficiency, acceptance angle etc.), as well 
as features specific to the target. This includes, in particular, the efficiency of contact 
with the target medium by the detector or intervening optical fibers that deliver and 
collect the optical signal. This is necessary because all model-based imaging schemes 
proposed for imaging in highly scattering media assume equivalency in detector 
efficiency for measured and computed values. Failure to account for such variables will 
introduce error whose magnitude and stability could vary considerably depending on the 
specifics of the measuring device and target medium. 

In principle, these uncertainties can be taken into account, although not without 
considerable effort for instrument design and added complexity of the forward modeling 
code. Generally speaking, such limitations are widely appreciated by many, both in this 
and other imaging communities. The goal in identifying practical schemes is to devise 
strategies that are mainly insensitive to such uncertainties. Often a desirable starting 
point is to employ schemes that provide useful information based on some type of 
relative measurement. Previously, we described a back projection scheme that evaluated 



28 



WO 01/20546 PCT/US00/25156 

relative detector data. RX. Barbour, H. Graber, R. Aronson, and J. Lubowski, "Model for 
3-D optical imaging of tissue " Int. Geosci and Remote Sensing Symp., (IGARSS), 2, 
1395-1399 (1990) (the disclosures of which are incorporated herein by reference). H.L. 
Graber, J. Chang, J. Lubowsky, R. Aronson and R.L. Barbour, "Near infrared absorption 
5 imaging in dense scattering media by steady-state diffusion tomography", in Proc. 

Photon migration and imaging in random media and tissues 1888, 372-386, (1993) (the 
disclosures of which are incorporated herein by reference). This formulation employed 
model-based imaging operators, but produced solutions lacking physical units. The lack 
of physical units (1) makes specific interpretation difficult, especially should multi- 
Q 10 wavelength measurements be considered, and (2) makes efforts to compute iterative 
00 updates difficult. However, this scheme however showed that in all cases tested, high 
M contrast images having excellent edge detection and object location were achievable. R.L. 
FU Barbour, H. Graber, R. Aronson, and J. Lubowski, "Model for 3-D optical imaging of 
J^j tissue," Int. Geosci and Remote Sensing Symp., (IGARSS), 2, 1395-1399 (1990) (the 
J 15 disclosures of which are incorporated herein by reference). H.L. Graber, J. Chang, J. 
pi Lubowsky, R. Aronson and R.L. Barbour, "Near infrared absorption imaging in dense 
scattering media by steady-state diffusion tomography", in Proc. Photon migration and 
imaging in random media and tissues", 1888, 372-386, (1993> (the disclosures of which 
are incorporated herein by reference). H.L. Graber, J. Chang and R.L. Barbour "Imaging 
20 of multiple targets in dense scattering media", in Proc. Experimental and numerical 

methods for solving ill-posed inverse problems: Medical and non-medical applications, 
SPIE, 2570, 219-234, (1995) (the disclosures of which are incorporated herein by 
reference). 
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In subsequent studies we have identified that the expression evaluated using this 
method has a functional form of: 



element number and N is the number of elements. It is apparent that this expression is 
similar to the following equation, which is equation (3) with a different form, 



While it is evident the two expressions are not equivalent, a more careful 
examination reveals that the different quantities on the left hand side of Eqs. (8) and (9) 
(i.e.,, the sum of weights and the reference intensity for the i th source-detector pair) are 
closely related. The influence that different forms of the data vector have on image 
recovery is discussed subsequently. 

Solution of the perturbation formulation requires specification of three input data 
sets. Two quantities, I r and W r , are typically computed from a specified reference 
medium, and the third quantity is the measured response, I. Results presented in FIGS. 6 
through 9 have explored the influence that errors in each of these quantities have on the 
quality of computed reconstructions for the original perturbation formulation (equation 
(1)). Most sensitive was the case where both quantities I r and W r are in error (c.f, test 
case 3 and FIG. 6). This suggests that the origin of excessive ill-conditioning can be 

M 

traced to one or both of the quantities I r and W r . When prior knowledge W r is assumed, 
errors in I r (test case 6) still produced highly unstable solutions (FIG. 8). On the other 
hand, when prior knowledge of the correct values of I r is assumed (test case 5), the 




(8) 



for each source-detector pair where, i is the source-detector pair number, j is the 



(9) 
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sensitivity to errors in W r is much less (FIG. 7), at least qualitatively. This strongly 
suggests that the principal cause of instability can be traced to errors in I r . It also 
indicates that should error free measurement data be available, significant instability in 
the solution domain would persist. The influence of systematic error in measurement was 
explored in test case 7. Results in FIG. 9 showed that under conditions where prior 
knowledge of I r and W r for the background is available, qualitatively accurate solutions 
could be obtained even in the presence of 900% error. In the absence of this prior 
knowledge, additional error due to this quantity further degraded image quality even in 
those cases where the selected reference medium differed only minimally from the 
background. 

While we have demonstrated that the described reconstruction procedure can 
significantly stabilize the computed reconstructions to errors in the reference medium, 
and that the principal origin of solution instability in its absence is the result of to errors 
in I r , it is useful to gain further insight as to why this should occur. An important 
difference between the two perturbation formulations is that equation (1) computes the 
difference between two exponentially attenuated quantities, while equation (3) performs a 
linear operation on an exponentially attenuated quantity. Because of the non-linear 
relationship between the medium coefficient values and surface detector responses, small 
errors in the former (the selected reference medium) can lead to large errors in the latter 
(the computed intensity or weight associated with the selected reference medium). 
Moreover, because the relationship is non-linear, such errors can be expected -to 
effectively distort the information content of the resultant data vector. The occurrence of 
such distortion is shown in FIGS. 16 and 17. FIG. 16 shows the angle dependence of the 
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relative detector response (i.e., the data vector) for a source bisecting the two inclusions 
computed using reference media corresponding to row 3 and column 3 of FIG. 5. 
Inspection of these plots clearly reveal a bimodal attenuation profile indicating the 
presence of two buried objects, a finding consistent with its actual structure. In contrast, 
this structure is almost completely absent from results derived using the original 
perturbation formulation (see FIG. 17) even though the variation in range of the reference 
medium is much less than that used in FIG. 16. This difference is most evident in results 
given in FIG. 1 8 that shows the amplitude of the frequency spectrum of the 
corresponding Fourier transforms. Comparison reveals that in the case where the selected 
reference medium has an error only in fz a of 0.02 cm" 1 , this produces an error of 
approximately 5 orders of magnitude in the amplitude of the frequency spectrum. This 
error grows to approximately eight orders of magnitude for the case where the selected 
reference medium has an error in fi s of only 5 cm' 1 . Given the magnitude of these errors, 
the extreme sensitivity of the reconstruction results obtained from the original 
perturbation equation is evident. 

When we examined the computed data vector for the formulation based on 
equation (7), we observed a pattern similar to that shown in FIG. 16, but with much 
higher amplitude values. It is this close relationship that we believe accounts for why 
previous reconstruction results derived using equation (8) also provided stable and 
qualitatively accurate maps, even though the solutions lacked features important for 
specific interpretation and iterative updates. R.L. Barbour, H. Graber, R. Aronson, and J. 
Lubowski, "Model for 3-D optical imaging of tissue," Int. Geosci and Remote Sensing 
Symp., (IGARSS), 2, 1395-1399 (1990) (the disclosures of which are incorporated herein 
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by reference). HX. Graber, J. Chang, J. Lubowsky, R. Aronson and RX. Barbour, "Near 
infrared absorption imaging in dense scattering media by steady-state diffusion 
tomography", in Proc. Photon migration and imaging in random media and tissues", 
1888, 372-386, (1993) (the disclosures of which are incorporated herein by reference). 
YX. Graber, J. Chang and RX. Barbour "Imaging of multiple targets in dense scattering 
media", in Proc. Experimental and numerical methods for solving ill-posed inverse 
problems: Medical and non-medical applications, SPIE, 2570, 219-234, (1995) (the 
disclosures of which are incorporated herein by reference). 

A further advantage of the current method compared to the previously described 
S ART-Type algorithm is that we are able to directly evaluate intensity difference values 
for which no mismatches exist between the computed intensity values, I r , and the 
Jacobian matrix, W r . YX. Graber, J. Chang and RX. Barbour "Imaging of multiple 
targets in dense scattering media", in Proc. Experimental and numerical methods for 
solving ill-posed inverse problems: Medical and non-medical applications, SPIE, 2570, 
219-234, (1995) (the disclosures of which are incorporated herein by reference). Not only 
does this reduce systematic errors, but produces solutions that can be updated by iterative 
methods. As reported here, error analysis studies have demonstrated that the current 
methodology produces remarkable stable solutions having excellent qualitative accuracy. 

An important goal shared by many investigators in the biomedical optics field is 
the capacity to accurately quantify variations in optical coefficients in tissue. One 
consequence of adopting the described reconstruction procedure is that the derived 
solution will not converge to the actual value even with error free measurement data and 
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non-linear updates. Instead, as is evident from results shown in FIG. 4 and 5, they will 
converge to a solution that is proportional to the true value throughout the image map. 

Examination of the recovered values indicates that the value of this 
proportionality coefficient depends strongly on the selected reference medium. 
Interestingly, though, as shown in FIG. 19, the ratio of the computed coefficient values is 
nearly constant over a broad range indicating that the relative error in the coefficients is 
itself a constant. In fact, for the examples studied, we find that (Sfia/dD) Reconstructed ss 
(dfiJdDf onginai over most of the range of reference values explored. In other examples 
(results not shown), we have explored this relationship for a variety of perturbation 
values with a similar range of reference media. In all cases, a constant error in the ratio of 
the derived coefficients was obtained. The value of the proportionality constant, 
however, varied depending on the magnitude and direction of the perturbation, but 
remained within a relatively small range. 

While it may be that case for some studies the restrictions imposed by the 
described method may prove limiting, it is anticipated that there are many practical 
situations where measurement of relative changes are nonetheless very useful. One 
example of special interest is in the evaluation of dynamic imaging data. 

Finally, from a mathematical perspective, the described methodology has the 
effect of desensitizing the solution of boundary value problems to specific features of the 
boundary. One practical consequence of this is that, unlike the standard perturbation 
formulation, the current formulation is less sensitive to a detailed knowledge of the 
external boundary of tissue, a quantity not easily measured. 
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While the instant invention has considered using of DC sources, it is expressly 
understood that the described methodology is extendable to other source conditions (e.g., 
time-harmonic and time-resolved measurements) and other inverse formulations (e.g., 
iterative gradient descent) and is extendable to other imaging modalities including, 
ultrasound, radioscintigraphic and impedance imaging. 

Although illustrative embodiments have been described herein in detail, those 
skilled in the art will appreciate that variations may be made without departing from the 
spirit and scope of this invention. Moreover, unless otherwise specifically stated, the 
terms and expressions used herein are terms of description and not terms of limitation, 
and are not intended to exclude any equivalents of the system and methods set forth in the 
following claims. 
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What is claimed is: 

1 . A method for imaging of the properties of a scattering target medium, 
comprising: 

generating a first vector of measured data and a second vector of measured data, 
the first vector of measured data being indicative of energy emerging from a target 
medium, the second vector of measured data being indicative of energy emerging from a 
target medium, the emerging energy substantially originating from at least one source 
directing the energy into the target medium; 

normalizing the first and second vectors of measured data; and 
solving a modified perturbation formulation of the radiation transport inverse 
problem for a relative change between a known property of a reference medium and the 
corresponding unknown property of a target medium, wherein the modified perturbation 
equation relates the normalized measured data and a vector of reference data for the 
known reference medium to the relative change in the property, the vector of reference 
data being indicative of energy emerging from the known reference medium. 

2. The method of claim 1 wherein the normalization of the first and second 
vectors of measured data comprises determining the difference between the first and 
second vectors of measured data relative to the second vector of measured data. 

3. The method of claim 1 wherein the modified perturbation equation has the 
following form: 
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where 8x is a vector of the relative changes between a known property of the reference 
medium and the corresponding unknown property of a target medium, for corresponding 
volume elements of the reference medium and the target medium, the volume elements 
being an imaginary grid of contiguous regions forming a representation of the target 
medium and reference medium, W\ is a weight matrix describing the influence that each 
of a plurality of volume elements of the reference medium has on energy emerging at a 
point on the reference medium, I r is the vector of reference data indicative of energy 
emerging from the reference medium, I is the first vector of measured data and I 0 is the 
second vector of measured data. 

4. The method of claim 1 wherein the normalization of the first and second 
sets of measured data comprise determining the natural logarithm of the quotient of the 
first set of measured data and the second set of measured data. 

5. The method of claim 1 wherein the modified perturbation equation has the 
following form: 
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where Sx is a vector of the relative changes between a known property of the 
reference medium and the corresponding unknown property of a target medium for 
corresponding volume elements of the reference medium and the target medium, the 
volume elements being an imaginary grid of contiguous, nonoverlapping regions forming 
a representation of the target medium and reference medium, W r is a weight matrix 
describing the influence that each of a plurality of volume elements of the reference 
medium has on energy emerging at a point on the reference medium, where I r is the 
vector of reference data indicative of energy emerging from the reference medium, I is 
the first vector of measured data and I 0 is the second vector of measured data. 

6. The method of claim 1 wherein the property is at least one of an 
absorption coefficient and a scattering coefficient. 

7. The method of claim 1 wherein the first vector of measured data and 
second vector of measured data are obtained from one target. 

8. The method of claim 1 wherein the first vector of measured data is 
obtained from a first target and the second vector of measured data is obtained from a 
second target. 

9. The method of claim 1 wherein the first vector of measured data is 
obtained at a first instant in time and the second vector of measured data is obtained at a 
second instant in time. 
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10. The method of claim 1 wherein the first vector of measured data is 
obtained at a first instant in time and the second vector of measured data is a time 
averaged mean of a plurality of measurements. 

11. A method of claim 1 further comprising generating an image representing 
the cross-sectional relative changes in the property. 

12. A system for imaging of the properties of a scattering target medium, 
comprising: 

means for generating a first vector of measured data and a second vector of 
measured data, the first vector of measured data being indicative of energy emerging 
from a target medium, the second vector of measured data being indicative of energy 
emerging from a target medium, the emerging energy substantially originating from at 
least one source directing the energy into the target medium; 

means for normalizing the first and second vectors of measured data; and 
means for solving a modified perturbation formulation of the radiation transport 
inverse problem for a relative change between a known property of a reference medium 
and the corresponding unknown property of a target medium, wherein the modified 
perturbation equation relates the normalized measured data and a vector of reference data 
for the known reference medium to the relative change in the property, the vector of 
reference data being indicative of energy emerging from the known reference medium. 
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13. The system of claim 12 wherein the normalization of the first and second 
vectors of measured data comprises determining the difference between the first and 
second vectors of measured data relative to the second vector of measured data. 

14. The system of claim 12 wherein the modified perturbation equation has 
the following form: 

W -Sx = 6\ , 
r r 

where 5x is a vector of the relative changes between a known property of the reference 
medium and the corresponding unknown property of a target medium, for corresponding 
volume elements of the reference medium and the target medium, the volume elements 
being an imaginary grid of contiguous regions forming a representation of the target 
medium and reference medium, W' r is a weight matrix describing the influence that each 
of a plurality of volume elements of the reference medium has on energy emerging at a 
point on the reference medium, I r is the vector of reference data indicative of energy 
emerging from the reference medium, I is the first vector of measured data and Io is the 
second vector of measured data. 

15. The system of claim 12 wherein the normalization of the first and second 
sets of measured data comprise determining the natural logarithm of the quotient of the 
first set of measured data and the second set of measured data. 
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16. The system of claim 12 wherein the modified perturbation equation has 
the following form: 

(w y 

K rhj , 

sv = w; 

where 8x is a vector of the relative changes between a known property of the 
reference medium and the corresponding unknown property of a target medium for 
corresponding volume elements of the reference medium and the target medium, the 
volume elements being an imaginary grid of contiguous, nonoverlapping regions forming 
a representation of the target medium and reference medium, W r is a weight matrix 
describing the influence that each of a plurality of volume elements of the reference 
medium has on energy emerging at a point on the reference medium, where I r is the 
vector of reference data indicative of energy emerging from the reference medium, I is 
the first vector of measured data and I 0 is the second vector of measured data. 
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(57) Abstract: A method for imaging of a scattering target medium using a modified perturbation formulation of the radiation trans- 
port equation wherein normalized measured values are used to recover a relative difference in absorption and/or scattering properties 
based on the normalized measured values with respect to a reference medium. The modified perturbation formulation provides en- 
hanced stability, reduces the sensitivity of solution to variations between the target and reference media, produces solutions having 
physical units and reduces the need for absolute detector calibration. Moreover, the modified perturbation equation lends itself to 
the detection and imaging of dynamic properties of the scattering medium. 
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OR PCT INTERNATIONAL APPLICATIONS (DESIGNATING THE U.S/> 

I hereby claim the benefit under Title 35, United States Code § 120 of any United States application(s) or under 
§ 365(c) of any PCT international application(s) designating the U.S. listed below. 



US/PCT Application Serial No. Filing Date Status (patented, pending, abandoned)/ 

U.S. application no. assigned (For PCT) 



US/PCT Application Serial No. Filing Date Status (patented, pending, abandoned)/ 

U.S. application no. assigned (For PCT) 



[ ] In this continuation-in-part application, insofar as the subject matter of any of the claims of this 
application is not disclosed in the above listed prior United States or PCT international application(s) in the manner 
provided by the first paragraph of Title 35, United States Code, § 1 12, 1 acknowledge the duty to disclose material 
information as defined in Title 37, Code of Federal Regulations, § 1.56(a) which occurred between the filing date of 
the prior application(s) and the national or PCT international filing date of this application. 

I hereby declare that all statements made herein of my own knowledge are true and that all statements made on 
information and belief are believed to be true; and further that these statements were made with the knowledge that 
willful false statements and the like so made are punishable by fine or Imprisonment, or both, under Section 1001 of 
Title 18 of the United States Code and that such willful false statements may jeopardize the validity of the 
application or any patent issued thereon. 
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I hereby appoint the following attorneys and/or agents with full power of substitution and revocation, to prosecute 
this application, to receive the patent, and to transact all business in the Patent and Trademark Office connected 
therewith: John A. Diaz (Reg. No. 19.550 ), John C. Vassil (Reg. No. 19,098). Alfred P. Ewert (Reg. No. 19,887 ), 
David H. Pfeffer (Reg. No. 19,825 ), Harry C. Marcus (Reg. No. 22,390), Robert E. Paulson (Reg. No. 21.046 ). 
Stephen R. Smith (Reg. No. 22,615), Kurt E. Richter (Reg. No. 24,052), J. Robert Dailey (Reg. No. 27,434 ). Eugene 
Moroz (Reg. No. 25,237 ), John F. Sweeney (Reg. No. 27.471 ). Arnold I. Rady (Reg. No. 26,601 ). Christopher A. 
Hughes (Reg. No. 26,914 ), William S. Feiler (Reg. No. 26.728 ). Joseph A. Calvaruso (Reg. No. 28.287 ). James W. 
Gould (Reg. No. 28,859), Richard C. Komson (Reg. No. 27,913 ), Israel Blum (Reg. No. 26,710) . Bartholomew 
Verdirame (Reg. No. 28,483) , Maria C.H. Lin (reg. No. 29.323) . Joseph A. DeGirolamo (Reg. No. 28,595 ), Michael 
P. Dougherty (Reg. No. 32,730) , Seth J. Atlas (Reg. No. 32,454 ), Andrew M. Riddles (Reg. No. 31.657 ). Bruce D. 
DeRenzi (Reg. No. 33,676 ), Michael M. Murray (Reg. No. 32.537 ). Mark J. Abate (Reg. No. 32.527 ). John T. 
3>^> Gallagher (Reg. No. 35,516), Steven F. Meyer (Reg. No. 35,613) , Kenneth H. Sormenfeld (Reg. No. 33.285 ). Tony 
V. Pezzano (Reg. No. 38.271) , Andrea L. Wayda (Reg. No. 43.979) and Walter G. Hanchuk Reg. No. ( 35,179 ) of 
Morgan & Finnegan, L.L.P. whose address is: 345 Park Avenue. New York, New York, 10154i and M ichael S. 
Marcus (Reg. NoTj3JLJ27) and John E. Hoel (Reg. No. 26*222) of Morgan & Finnegan, L.L.P., whose address is 
1775 Eye Street, Suite 400, Washington, D.C. 20006. 

[ ] I hereby authorize the U.S. attorneys and/or agents named hereinabove to accept and follow instructions 





regarding this application without direct communication between the U.S. attorneys and/or agents and me. 
In the event of a change in the person(s) from whom instructions may be taken I will so notify the U.S. 
attorneys and/or agents hereinabove. 



I _ Full name of sole or first inventor Randall L. BARBOUR 



date 

Residence 15 Cherry Lane, Glen Head, New York 1 1546 Sj V 




Citizenship United States 



Post Office Address 



same 



ru 



Full name of second joint inventor, if any 



Inventor's signature* 



date 



Residence 



Citizenship 



Post Office Address 
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Full name of third joint inventor, if any 

Inventor's signature* 

date 

Residence 

Citizenship 

Post Office Address 



[ ] ATTACHED IS/ARE ADDED PAGE(S) TO COMBINED DECLARATION AND POWER OF 
ATTORNEY FORM FOR SIGNATURE BY FOURTH AND SUBSEQUENT INVENTORS 

* Before signing this declaration, each person signing must: 

1 . Review the declaration and verify the correctness of all information therein; and 

2. Review the specification and the claims, including any amendments made to the claims. 

After the declaration is signed, the specification and claims are not to be altered. 
To the inventor(s): 

The following are cited in or pertinent to the declaration attached to the accompanying application: 

Title 37, Code of Federal Regulation, §1.56 

Duty to disclose information material to patentability. 

(a) A patent by its very nature is affect with a public interest. The public interest is best served, and 
the most effective patent examination occurs when, at the time an application is being examined, the Office 
is aware of and evaluates the teachings of all information material to patentability. Each individual 
associated with the filing and prosecution of a patent application has a duty of candor and good faith in 
dealing with the Office, which includes a duty to disclose to the Office all information known to that 
individual to be material to patentability as defined in this section. The duty to disclose information exists 
with respect to each pending claim until the claim is canceled or withdrawn from consideration, or the 
application becomes abandoned. Information material to the patentability of a claim that is canceled or 
withdrawn from consideration need not be submitted if the information is not material to the patentability 
of any claim remaining under consideration in the application. There is no duty to submit information 
which is not material to the patentability of any existing claim. The duty to disclose all information known 
to be material to patentability is deemed to be satisfied if all information known to be material to 
patentability of any claim issued in patent was cited by the Office or submitted to the Office in the manner 
prescribed by §§ 1 .97(b)-(d) and 1 .98. However, no patent will be granted on an application in connection 
with which fraud on the Office was practiced or attempted or the duty of disclosure was violated through 
bad faith or intentional misconduct. The Office encourages applicants to carefully examine: 

(1) prior art cited in search reports of a foreign patent office in a counterpart application, and 
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(2) the closest information over which individuals associated with the filing or prosecution of 
a patent application believe any pending claim patentably defines, to make sure that any 
material information contained therein is disclosed to the Office. 



Title 35. U.S. Code $ 101 
Inventions patentable 

Whoever invents or discovers any new and useful process, machine, manufacture, or composition of 
matter, or any new and useful improvement thereof, may obtain a patent therefor, subject to the conditions and 
requirements of this title. 



Title 35 U.S. Code § 102 

Conditions for patentability; novelty and loss of right to patent 
A person shall be entitled to a patent unless - 

(a) the invention was known or used by others in this country, or patented or described in a printed 
publication in this or a foreign country, before the invention thereof by the applicant for patent, 

(b) the invention was patented or described in a printed publication in this or foreign country or in 
public use or on sale in this country, more than one year prior to the date of application for patent in the United 
States, or 

© (c) he has abandoned the invention, or 

jM. 

CI (d) the invention was first patented or caused to be patented, or was the subject of an inventor's 

|1 1 certificate, by the applicant or his legal representatives or assigns in a foreign country prior to the date of the 
fi . application for patent in this country on an application for patent or inventor's certificate field more than twelve 
q months before the filing of the application in the United States, or 

jj (e) the invention was described in a patent granted on an application for patent by another filed in the 

United States before the invention thereof by the applicant for patent, or on an international application by another 
3J wno nas fulfilled the requirements of paragraphs (1), (2), and (4) of section 371(c) of this title before the invention 
5f thereof by the applicant for patent, or 

in 

(f) he did not himself invent the subject matter sought to be patented, or 

(g) before the applicant's invention thereof the invention was made in this country by another had not 
abandoned, suppressed, or concealed it. In determining priority of invention there shall be considered not only the 
respective dates of conception and reduction to practice of the invention, but also the reasonable diligence of one 
who was first to conceive and last to reduce to practice, from a time prior to conception by the other . . . 



Title 35, U.S. Code $ 103 

Conditions for patentability; non-obvious subject matter 

A patent may not be obtained though the invention is not identically disclosed or described as set forth in 
section 102 of this title, if the differences between the subject matter sought to be patented and the prior art are such 
that the subject matter as a whole would have been obvious at the time the invention was made to a person having 
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ordinary skill in the art to which said matter pertains. Patentability shall not be negatived by the manner in which 
the invention was made. 

Subject matter developed by another person, which qualifies as prior art only under subsection (f) or (g) of 
section 102 of this title, shall not preclude patentability under this section where the subject matter and the claimed 
invention were, at the time the invention was made, owned by the same person or subject to an obligation of 
assignment to the same person. 



Title 35, U.S. Code § 1 12 Onparf) 
Specification 

The specification shall contain a written description of the invention, and of the manner and process of 
making and using it, in such full, clear, concise and exact terms also enable any person skilled in the art to which it 
pertains, or with which it is mostly nearly connected, to make and use the same, and shall set forth the best mode 
contemplated by the inventor of carrying out his invention. 



E3 



ni 



P£5 
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Title 35, U.S. Code $ 119 

Benefit of earlier filing date in foreign country; right of priority 

An application for patent for an invention filed in this country by any person who has, or whose legal 
representatives or assigns have, previously regularly filed an application for a patent for the same invention in a 
foreign country which affords similar privileges in the case of applications filed in the United States or to citizens of 
the United States, shall have the same effect as the same application would have if filed in this country on the date 
on which the application for patent for the same invention was first filed in such foreign country, if the application in 
this country is filed within twelve months from the earliest date on which such foreign application was filed; but no 
patent shall be granted on any application for patent for an invention which had been patented or described in a 
printed publication in any country more than one year before the date of he actual filing of the application in this 
country, or which had been in public use or on sale in this country more than one year prior to such filing. 

Title 35, U.S . Code § 120 

Benefit or earlier filing date in the United States 

An application for patent for an invention disclosed in the manner provided by the first paragraph of section 
1 12 of this title in an application previously filed in the United States, or as provided by section 363 of this title, 
which is filed by an inventor or inventors named in the previously filed application shall have the same effect, as to 
such invention, as though filed on the date of the prior application, if filed before the patenting or abandonment of or 
termination of proceedings on the first application or an application similarly entitled to the benefit of the filing date 
of the first application and if it contains or is amended to contain a specific reference to the earlier filed application. 

Please read carefully before signing the Declaration attached to the accompanying Application. 

If you have any questions, please contact Morgan & Finnegan, L.L.P. 



FORMrCOMB-DEC.NY 
Rev. 10/00 
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